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Abstract. A third post-Newtonian (3 PN) equation of motion for two spherical 
compact stars in a harmonic coordinate has been derived based on the surface 
integral approach and the strong field point particle limit. The strong field point 
particle limit enables us to incorporate a notion of a self-gravitating regular star 
into general relativity. The resulting 3 PN equation of motion is Lorentz invariant, 
unambiguous, and conserves an energy of the binary orbital motion. 
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1. Introduction 

Recently, high order post-Newtonian iteration for inspiralling compact binaries 
has attracted a renewed attention in the context of the worldwide gravitational 
waves detection efforts ^ |21 El E]- Detectability and qualities of measurements of 
gravitational waves emitted by such binaries depend on our theoretical knowledge 
about waveforms and hence partly about dynamics of binaries [S] . 

The currently available equation of motion and the Hamiltonian that govern 
dynamics of (comparable masses) binaries in the post-Newtonian approximation of 
general relativity are the 3.5 PN equations of motion in a harmonic coordinate [HIEHS] 
and the 3.5 PN Arnowitt-Deser-Misner (ADM) Hamiltonian in the ADM transverse- 
traceless (ADMTT) gauge j^l^l^J. For reviews on post-Newtonian iteration, we 
refer readers to, e.g., [T^j . 

In the inspiralling phase of spherical compact stars binaries, a point particle 
description is certainly useful and suitable to study their dynamics. A simple 
realization of this description may be a use of Dirac delta distributions, however, it 
causes ill-defined divergences due to the non-linearity of general relativity. The above 
mentioned works (except for [Sj) have used Dirac delta distributions and resorted 
to regularization procedures such as Hadamard's partie finie. Then a problem has 
occurred. At the 3 PN order, there exists one arbitrary undetermined parameter 
which reduces the predictability of the theoretical waveforms. Jaranowski and Schiifer 
[1UI IT3*] who first achieved the 3 PN iteration found two undetermined coefficients; 
Wk and tu s in their Hamiltonian (later they fixed Wk by imposing Poincare invariance 
on their Hamiltonian ^1]). Blanchet and Faye E] found one coefficient, called 
A, in their equation of motion. They then showed [151 ITfi| that the results of these 
two independent studies are consistent with each other when the following relation 
holds; lu s = — 11A/3 — 1987/840. It was suggested that unsatisfactory features of the 
regularizations they have used associated with their use of Dirac delta distributions 
were possible origins of the parameters, as summarized in |11| . 

In fact, Damour, Jaranowski, and Schafer succeeded in deriving the complete 
3 PN ADM Hamiltonian in the ADMTT gauge and fixed the parameter as ui s = 
using dimensional regularization [llj . However, applicability of mathematical 
regularizations to the current problem is not a trivial issue, but an assumption to 
be verified, or at least supported by convincing arguments. At the 2.5 PN order, 
there is a physical argument of the "Dominant Schwarzschild Condition" [17| which 
supports a use of Dirac delta distributions and a certain sort of regularizations. There 
is no such argument at the 3 PN order. Thus, it seems crucial to achieve unambiguous 
3 PN iteration without introducing singular sources and support the result of |11| . 

The author, Futamase and Asada have recently developed yet another method 
of derivation of a post-Newtonian equation of motion in a harmonic coordinate 
appropriate to inspiralling compact binaries [181 I19j which is based on the strong 
field point particle limit 20 and the surface integral approach [21] . In this paper, due 
to the lack of space, we briefly review our method and present few of our 3 PN results. 
Full explanations on our method and results will be reported elsewhere [221 EH ■ 

2. Strong field point particle limit 

We shall write the post-Newtonian expansion parameter explicitly as e, which 
represents smallness of the orbital velocity, O(e), and weakness of the interstars field, 
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0(e 2 ). Futamase [20] advocated the strong field point particle limit where the radius 
of the star shrinks at the same rate as the mass in the post-Newtonian limit (e — > 0). 
Thus, the typical scale of gravitational field inside the star (~ the mass over the radius 
of the star) remains finite (and can be strong) in e zero limit. This is in contrast with 
the usual post-Newtonian expansion for an extended star, where the matter density 
is assumed to scale as e 2 and we apply the post-Newtonian approximation to the field 
even inside the stars. 

3. Field equations 

Using the post-Newtonian approximation, we iteratively solve the Einstein equations 
under the harmonic gauge condition that is written as h^ v „ = where W = 
V ltu ~ V~99^ L ' i V 1 *" = diag(— e 2 , 1, 1, 1) is the flat metric, and g is the determinant 
of the metric g^ v . Then harmonically relaxed Einstein equations can be casted 
in integral forms as /^(r, x 4 ) = 4 $ c <ftyk tlv {T — e\x — y\,y k ;e)/\x — y\, where 
A MI/ = + x MI/a/3 ,a/3: 8 M1/ is the matter plus the gravitational field stress energy 
pseudo-tensor (8 piy = —g{T^ v + t^ L ) with t^ L being the Landau-Lifshitz pseudo 
tensor) and X^" ,ap arises since we use the wave operator of the flat spacetime instead 
of that of the curved spacetime. r is the time coordinate. The light cone C emanating 
from the field point (r, x l ) is split into two zones; the near zone N which encloses 
the two stars completely and the far zone F which is the outside of N and where 
the retardation effects of the field are manifest. To treat the retarded integrals over 
F we have adopted the Direct Integration of the Relaxed Einstein equations method 
developed by Will and collaborators jS], however, it is well-known jSJ |^] that the 
integrals over F do not contribute to an equation of motion up to the 3.5 PN order 
and hence we devote ourselves to the integrals over TV. 

The integral region N then is split into three regions; two spheres, called the 
body zone Ba (A = 1,2 labels each star), each of which surrounds each star with no 
intersection with the other nor the star, and elsewhere N/B where B = B\ U B^. We 
make the radius of Ba shrink proportionally to e to define the multipole moments and 
to derive an equation of motion for compact stars. That is, Ba = {x l \\x — z*a(t)\ < 
cRa} where cRa is arbitrary but smaller than the orbital separation r\2 and larger 
than the radius of the star for any e. z\(t) represents the location of the star A in 
the point particle limit. The integrals over each Ba are evaluated using the multipole 
expansion of the star. For example, the (8 part of the) mass multipole moments of 

the star A are defined as ify = e 2 f B d 3 aAQ TT ct^ L , where or A = e~ 2 (x l — z\{t)) is 
a scaled coordinate in which the star does not shrink to the point z A (r). The capital 

index denotes a set of collective indices, I\ = ix ■ • ■ ii and crj = cr^a^f ■ ■ ■ a a- 

The remaining Poisson-type integrals over N/B, J N , B d^y\x— y\ n ~ 1 f(r, y) (n > : 
integer, f(r, x) being a certain combination of A M!y (r, x)), are evaluated with the help 
of particular solutions of Poisson equations. For n — case, using A.g(x) = f(x) (for 
brevity, we do not write r dependence in / nor g here), we have 

/ ^M = -Ang(g) + i dS k 

J N/B F - V\ Jd(N/B) 

8nce explicit expressions of the necessary particular solutions are in hand, it is 
straightforward to evaluate the Poisson-type integrals. 



1 dg{y) 
\x-y\ dy k 
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4. Mass-energy relation and momentum-velocity relation 

In our formalism, we do not assume a specific expression of the matter's stress energy 
tensor. Thus, we need to express the operationally defined multipole moments, 
namely, the four momentum P A q = e 2 J B d 3 <XA@ >J ' T of the star in terms of the 
mass vtla, the velocity va, and r%2. As for P A q, we use its evolution equation which 
is derived by the local energy momentum conservation 0*"" v — as dP A& /dr = 
—£~ 4 § dBA dSk{—g)t k L L + e~ 4: v A § dBA dSk{—g)tYL- At the lowest order, we obtain 
dP AQ /dr = 0(e 2 ). The evolution equation of P\& is then solved functionally into 
a form P Ae = m^fl + 0(e 2 )] where the mass of the star A, tua, is defined as the 
integration constant by vtla = lim e ^o Pas ■ We note that although uia is defined by 
e — * limit of P A& , the gravitational field energy inside the star is taken into account 
in the definition of vtla since Pas is defined as a volume integral whose integrand 
includes the field energy {—g)t T jJ L which does not vanish in the e zero limit because 
of the scaling of the strong field point particle limit, (the mass over the radius of the 
star is finite inside the star independently of e). Following the Newtonian dynamics, 
the momentum- velocity relation on the other hand is derived by the time derivative of 
the dipole moment D Ae = e 2 J B d 3 aA® TT a A as P A@ = Paq v a + Qa& + e 2 dD Ae /dr 
where Q* Ae = e" 4 /^ dS k ((-g)r L k L -v k A (-g)t L l)y A = -e e d(m AO a A /6)/dT ( a A is the 
Newtonian acceleration) arises due to the non-compactness of (—g)t^ L . Q A q appears 
starting from the 3 PN order and it affects the 3 PN equation of motion. Finally, we 
note that as in the Newtonian case, we can choose the value of the dipole moment 
D Ae freely to define the representative point of the star A, z a (t). For instance, z a (t) 
corresponding to D l Ae — may be called the center of mass of the star A. However, 
we found that it is convenient to adopt D Ae = e 4 (l/6 — 22 \n(ri 2 / e / RA))m A oa A at 
the 3 PN order [221 Choosing a certain value for D Ae corresponds to choosing the 
corresponding coordinate under the harmonic gauge condition. Our particular choice 
gauges \neRA dependence away from the 3 PN acceleration. A different choice leaves 
In sRa dependences in an acceleration, however, these are mere gauge terms and do 
not appear when we are concerned with coordinate invariant quantities. In fact, a 
different choice P>\q — e rn? A Qa A /& gives the same invariant binary orbital energy of 
circular orbital motion in the center of mass frame as E(x) given below |23|. 

5. Equation of motion 

There may be three methods being studied in derivation of a high order post- 
Newtonian equation of motion. The first one is the volume integral approach, 
which may take the following form in the Newtonian dynamics: mAdv A /dr = 
J d 3 ypA(y)4>(y)' 1 ■ (pa denotes the density of the star A. <p is the Newtonian potential.) 
In this method, one can track how the stars' internal structures dependent terms cancel 
out (or remain) and disappear (or appear) in the final result; if the resulting equation 
of motion depends on, say, the density profile of (almost stationary) spherically 
symmetric compact stars, it suggests violation of the equivalence principle to some 
extent. This method is being used by Pati and Will ; 8 ; . 

The second method assumes that a regularized action describes dynamics of 
binaries. Blanchet and Faye |5J [7] have used this method with the generalized 
Hadamard's partie finie regularization but obtained one arbitrary parameter A in 
their 3 PN equation of motion. It is shown jSJ [7] that a regularized geodesies can be 
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derived from a certain action. This method is attractive, besides its mathematical 
beauties, because if the resulting equation coincides with the equations derived by 
the first (and the third below) method, then it suggests a notion of "geodesies in a 
dynamical spacetime" (with a certain class of regularizations supplemented). 

The third method, which we have adopted, is the surface integral method where 
we compute the gravitational momentum flux t 13 going through a sphere surrounding 
the star concerned. In the Newtonian case, we may have mAdv l A / dr — — J dSjt lJ ; 
t % * = l/(47r)(</>' l J — <5 y 0' fc : fc/2). In this method, possible information of the stars' 
internal structures are coded in, say, the multipole moments of the stars and hence 
it affects the equation of motion through t %3 . Based on the local energy momentum 
conservation, we obtain in our formalism a simple generalization of the surface integral 
approach to general relativity as 



m A ^ = - i dS k (-g)1% L +v k A i dS k (-g)t T L 



-iaB A JOB A 



LL 



+ ( m ,-p; e )$ + ^-%-^k. (1) 

6. Gravitational field and equations of motion up to the 3 PN order 

Up to the 2.5 PN order, it is possible to derive the gravitational field explicitly by the 
method described above. Once we have the field in closed form, it is straightforward 
to evaluate the surface integrals in Eq. Q). However, it does not seem possible to 
find explicitly all the particular solutions for the Poisson-type integrals over N/B 
necessary to evaluate the 3 PN gravitational field, even in the neighborhood of the 
star. We thus partly abandon derivation of the field and took an alternative method 
22 . We change the order of integration; for the integrands for which we could not find 
the particular solutions of Poisson equations, we first evaluate the surface integrals in 
Eq. Q before evaluating the Poisson-type integrals and then evaluate the remaining 
volume integrals. With this method, we obtain a 3 PN equation of motion. Because 
of the enormous length of the final result, we here present only the 3 PN relative 
acceleration in the case of the circular orbit and in the center of mass frame, which is 
an appropriate equation to inspiralling binaries; dV /dr = — £l 2 r\ 2 + e 5 2.5PN^4' where 
V 1 = v\ — v\ is the relative velocity and 2.5PN^4 l is the relative acceleration at the 2.5 
PN order (the radiation reaction term). The 3 PN orbital angular frequency £1 is, 



m 2 Q 2 = 7 3 



1 + e 2 7 (-3 + v) + e V \6+^-v + v 2 ) r'VM-KI 



2375 . 41tt 2 1 19 2 3 



24 64 



-V + V 



G(e 7 ), (2) 



where m = mi + m,2, v = mimj/m 2 , 7 = mjvxi- We here note that it is not allowed 
to fix the A parameter by comparing Eq. @ with the corresponding result of Blanchet 
and Faye since the harmonic gauge condition both groups have used does not fix 
a coordinate completely |231 12*2] . 

Although we can not show it here due to the lack of space, our equation of motion 
respects the Lorentz invariance. And more importantly, our 3 PN equation of motion 
is unambiguous and admits a conserved energy. In the case of the circular orbit 
(neglecting the radiation reaction force) and in the center of mass frame, the invariant 
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orbital energy E becomes 

. . mi 

E{X) = ; 



,3 1 \ 4 / 27 19 1 M o 
! + ({--- —v \x + e 4 -— + — v - — v 2 x 2 



4 12 J V 8 8 24 



, ; , 675 f 34445 205^1 155 2 JS5_ 3 , .. 
64 + i 576 96~~ 1 ~~96~^ "5184^ 



where x = (mil) 2 / 3 . We have thus derived a 3 PN equation of motion free from an 
arbitrary parameter. Since we introduced Ra in the definition of Ba, it may seem 
that Ra appears in the field and also in the equation of motion. Ba is introduced to 
separate the region where the presence of the star may make the field strong from well 
outside of the stars where the post-Newtonian approximation to the field is applicable. 
The contribution to the field from Ba is evaluated with the multipole expansion. The 
equation of motion is derived via the surface integrals over OB a- Now, as for Ra 
dependent terms in the field, since we have used the non-singular matter sources, the 
field is smooth near 8Ba and possible Ra dependences in the body zone contributions 
to the field cancel out the corresponding Ra dependences in the N/B contributions. 
As for an equation of motion, we have proved in |19| that an equation of motion 
derived from Eq. does not depend on Ra- See f° r more details. 

Comparing our E(x) with that in we determine the coefficient undetermined 
in the Blanchet and Faye 3 PN equation of motion as A = —1987/3080. This value 
of A is consistent with the result of Damour, Jaranowski, and Schafer [11| . Thus, 
our result (indirectly) validates their use of dimensional regularization in the ADM 
Hamiltonian approach in the ADMTT gauge. 

In the paper |22l I23j . we devote ourselves to two spherical compact binaries, 
however, our formalism can be extended to include higher order multipole moments 
of the stars ^Hj- We have used the strong field point particle limit and the surface 
integral approach to treat strongly self-gravitating stars. Therefore our 3 PN equation 
of motion is applicable to regular compact stars with strong self-gravity. As is clear 
from Eq. J5J), our equation of motion depends only on tua, and this fact supports the 
strong equivalence principle up to the 3 PN order. 

Finally, as mentioned above, we could not derive the complete 3 PN gravitational 
field in closed forms. This fact makes it difficult to proceed to the 4 PN iteration. The 
enormous numbers of terms in the intermediate computations, of order of 10 5 , appear 
even at the 3 PN order. Brute force post-Newtonian iteration at the 4 PN order may 
need efficient computations techniques. It may be possible to adopt the strong field 
point particle limit in the currently available waveform computation formalism and 
fix the undetermined parameters in the 3.5 PN accurate waveforms. 
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